Consignes

Complétez ce document en remplissant les chunks vides pour écrire le code qui vous a permis de répondre à la question. Les réponses attendant un résultat chiffré ou une explication devront être insérés entre le balises html code. Par exemple pour répondre à la question suivante :

La bioinfo c'est : <code>MERVEILLEUX</code>.

N’hésitez pas à commenter votre code, enrichier le rapport en y insérant des résultats ou des graphiques/images pour expliquer votre démarche. N’oubliez pas les bonnes pratiques pour une recherche reproductible ! Nous souhaitons à minima que l’analyse soit reproductible sur le cluster de l’IFB.

Introduction

Vous allez travailler sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

# Creer son espace de travail
## Faire un répetoire pour chaque type de données
mkdir -p ~/M4-5-Eval/FASTQ  
mkdir -p ~/M4-5-Eval/CLEANING
mkdir -p ~/M4-5-Eval/MAPPING
mkdir -p ~/M4-5-Eval/QC
mkdir -p ~/M4-5-Eval/REF

# Vérification de l'espace de travail 
cd ~/M4-5-Eval
tree ~/M4-5-Eval

Téléchargement des données brutes

Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]

# Il faut en premier lieu loader le module sra-tools puis spécifier le répertoire de destination est FASTQ
module load  sra-tools
fasterq-dump -h
srun --cpus-per-task=6 fasterq-dump --split-files -p SRR10390685 --outdir FASTQ

# on zip les données 
cd ~/M4-5-Eval/FASTQ
srun gzip *.fastq

Combien de reads sont présents dans les fichiers R1 et R2 ?

#On sait que la synthaxe FASTQ est faite de telle sorte que 1 seul read s'écrit sur 4 lignes (Info avec numéro de machine, la séquence, le + qui est un séparateur et sur la dernière ligne le code de qualité du read). On compte donc le nombre de lignes présentes dans le fichier fastq et on le divise par 4. 

echo $(zcat SRR10390685_1.fastq.gz|wc -l)/4|bc
echo $(zcat SRR10390685_2.fastq.gz|wc -l)/4|bc

# on peut également utiliser le code suivant avec la commande seqkit (module qu'il faut également loader), cela permet de vérifier si les résultats sont les mêmes 
module load seqkit
srun seqkit stats --threads 1 *.fastq.gz

#file                    format  type   num_seqs        sum_len  min_len  avg_len  max_len
#SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
#SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151

Les fichiers FASTQ contiennent 7066055 reads.

Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

# se mettre dans le répertoire de travail REF : 

cd ~/M4-5-Eval/REF
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz

Quelle est la taille de ce génome ?


srun seqkit stats --threads 1 GCF_000009045.1_ASM904v1_genomic.fna 

La taille de ce génome est de 4,215,606 paires de bases.

Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz

Combien de gènes sont connus pour ce génome ?

# décompression du fichier gff.gz pour faire un awk dessus
gunzip GCF_000009045.1_ASM904v1_genomic.gff.gz

# awk pour sélectionner dans la colonne 3 ceux qui ont un identifiant 'gene' et avec wc -l je demande l'affichage du nombre de lignes. 
awk '$3=="gene"' GCF_000009045.1_ASM904v1_genomic.gff | wc -l

4448 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit


module load fastqc
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8

La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

  • Oui
  • Non

car On remarque que les valeurs de score qualité sont > 30 (> Q30). comme le montre le rapport de controle qualité Fastq et plus précisément le graphique “FastQ mean quality Scores” comme le montre

Lien vers le [rapport MulitQC] (https://fabienne-nr.github.io/Evaluation_M4M5/multiqc_report.html)

Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

  • Oui
  • Non

car les reads sont de différentes longueurs comme montrés sur le graphique Sequence Length Distribution

Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?

# Utilisation de la formule
# C = N*L/G ; avec N = Nombre de lectures total (ici 7066055 pour chaque reads donc il faut multiplier par 2), L = Longueur des lectures et G = Taille du génome 
# ((7066055 * 2) * 151) / 4215606 = 506.2021  

La profondeur de séquençage est de : 506 X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.

# Changer de repertoire de travail
cd ..

# charger le module fastp
module load fastp

# Récupérer les fastq.gz (en entrée --int1 et -- int2), spécifier la sortie dans le dossier CLEANING avec les options --out1 et --out2 et indiquer le fichier fastp.html et fastp.json dans le repertoir CLEANING. 

srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail --json CLEANING/fastp.json
Parametre Valeur Explication
cut_mean_quality 30 on accepte une erreur sur 1000 et on a donc une précision de base de 99.9%
cut_window_size 8 Spécifie le nombre de bases pour faire la moyenne de la qualité
Length_required 100 on ne garde que les reads de 100 pb minimum

Ces paramètres ont permis de conserver 13554096 reads pairés, soit une perte de 4,09% des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].

# changement de repertoire de travail 
cd ~/M4-5-Eval/MAPPING/

# chargement du module bwa
module load bwa

# Indexer la séquence de référence qui se trouve dans le repertoire REF 
srun bwa index ~/M4-5-Eval/REF/GCF_000009045.1_ASM904v1_genomic.fna.gz

# Faire un bwa mem
srun --cpus-per-task=32 bwa mem ~/M4-5-Eval/REF/GCF_000009045.1_ASM904v1_genomic.fna.gz ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685.sam

# Transformation du fichier sam en bam via samtools
module load samtools
srun --cpus-per-task=8 samtools view --threads 8 SRR10390685.sam -b > SRR10390685.bam
srun samtools sort SRR10390685.bam -o SRR10390685.sort.bam
srun samtools index SRR10390685.sort.bam

# Statistiques 
srun samtools idxstats SRR10390685.sort.bam > SRR10390685.sort.bam.idxstats
srun samtools flagstat SRR10390685.sort.bam > SRR10390685.sort.bam.flagstat

# multiQC
multiqc .d .o

Combien de reads ne sont pas mappés ?

cat ../samtools flagstat SRR10390685.sort.bam > SRR10390685.sort.bam.flagstat


#13571369 + 0 in total (QC-passed reads + QC-failed reads)
#0 + 0 secondary
#17273 + 0 supplementary
#0 + 0 duplicates
#12826829 + 0 mapped (94.51% : N/A)
#13554096 + 0 paired in sequencing
#6777048 + 0 read1
#6777048 + 0 read2
#12746370 + 0 properly paired (94.04% : N/A)
#12769290 + 0 with itself and mate mapped
#40266 + 0 singletons (0.30% : N/A)
#0 + 0 with mate mapped to a different chr
#0 + 0 with mate mapped to a different chr (mapQ>=5)

## Le nombre de reads non mappés 
# unmapped = "total" - "mapped" - "singletons"

echo "13571369-12826829-40266" | bc

704274 reads ne sont pas mappés.

Croisement de données

Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:

# chargement du module 
module load bedtools

# récupération du gène d'intérêt 
grep trmNF ../REF/GCF_000009045.1_ASM904v1_genomic.gff | awk '$3=="gene"' > trmNF_gene.gff3


# bedtools intersect avec -f à 0.5
srun bedtools intersect -a SRR10390685.sort.bam -b trmNF_gene.gff3 -sorted -f 0.5 > SRR10390685_trmNF.bam 


# Tri et indexage des reads
srun samtools sort SRR10390685_trmNF.bam -o SRR10390685_trmNF.sort.bam
srun samtools index SRR10390685_trmNF.sort.bam
srun samtools idxstats SRR10390685_trmNF.sort.bam > SRR10390685_trmNF.sort.bam.idxstats

# Calcul du nombre de read chevauchant 
cat SRR10390685_trmNF.sort.bam.idxstats 

#NC_000964.3    4215606 2801    0
#@* 0   0   0

2801 reads chevauchent le gène d’intérêt.

Visualisation

Utilisez IGV [7] sous sa version en ligne pour visualiser les alignements sur le gène. Faites une capture d’écran du gène entier.

Visualisation

References

1. toolkit NS. NCBI SRA toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.
4. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997. 2013.
5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
7. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (IGV): High-performance genomics data visualization and exploration. Briefings in bioinformatics. 2013;14:178–92.
 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France